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Introduction 

Movements of the celestial bodies in our solar system inspired Isaac Newton to work out his profound 
laws of gravitation and motion; with one or two notable exceptions, all of those objects move as Newton 
said they would. But normally harmonious orbital motion is accompanied by the risk of collision, which 
can be cataclysmic. The Earth’s moon is thought to have been produced by such an event, and we 
recently witnessed magnificent bombardments of Jupiter by several pieces of what was once Comet 
Shoemaker-Levy 9. Other comets or asteroids may have met the Earth with such violence that dinosaurs 
and other forms of life became extinct; it is this possibility that causes us to ask how the human species 
might avoid a similar catastrophe, and the answer requires a thorough understanding of orbital motion. 

The two red square flags with black square centers displayed in figure 1 are internationally recognized 
as a warning of an impending hurricane. Mariners and coastal residents who know the meaning of this 
symbol and the signs evident in the sky and ocean can act in advance to try to protect lives and property; 
someone who is unfamiliar with the warning signs or chooses to ignore them is in much greater jeopardy. 
Although collisions between Earth and large comets or asteroids occur much less frequently than landfall 
of a hurricane, it is imperative that we learn to identify the harbingers of such collisions by careful 
examination of an object’s path. 

An accurate determination of the orbit of a comet or asteroid is necessary in order to know if, when, 
and where on the Earth’s surface a collision will occur. Generally speaking, the longer the warning time, 
the better the chance of being able to plan and execute action to prevent a collision. The more accurate 
the determination of an orbit, the less likely such action will be wasted effort or, what is worse, an effort 
that increases rather than decreases the probability of a collision. Conditions necessary for a collision to 
occur are discussed, and warning times for long-period comets and near-Earth asteroids are presented. 

Orbit determination is the process of using a collection of measurements obtained by observation to 
calculate a set of orbital elements, six quantities that give (either implicitly or explicitly) the position and 


^Chapter nomenclature available in chapter notes, p. 217. 
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Figure 1. Flurricane Warning Flag. 


velocity of an object at a particular instant of time. The classical methods of Laplace, Gauss, and Olbers 
utilize six angular measurements, obtained two at a time from each of three optical observations, to make 
a preliminary determination of orbital elements. Additional measurements can be employed together with 
the iterative method of least squares to improve accuracy of the elements. With all these methods, the 
quality of the result is affected by errors in the measurements and by the spatial and temporal spacing of 
observations which are in turn related to the geometry resulting from the orbits of the object and the 
observatory. The method of least squares makes it possible to take advantage of multiple measurements 
made from a single observatory or from two or more observatories spread throughout the solar system; 
the accuracy of orbital elements so obtained is thus a function of the number and spacing of observations, 
and the number and placement of observatories. 

A study of the aforementioned factors and their effects on orbit determination is undertaken in order to 
identify trends and make a preliminary determination as to the number, placement, and resolution of 
instalments required to form an effective system for determining orbits of comets and asteroids. In addi- 
tion, an idea can be obtained of the number and timing of observations that yield the best results, leading 
to a strategy for using the system to perform observations. In order to keep this stage of the analysis from 
becoming unnecessarily complicated it is assumed that measurements can be made at will, but of course 
this cannot always be the case. Orbit determination of long-period comets using a batch filter and sequen- 
tial filters is examined, and a batch filter is applied to the study of near-Earth asteroids. 

Collisions 

This analysis is facilitated by the assumption that a collision occurs when the heliocentric distance of a 
comet or asteroid is identical with that of the Earth; in the case where the object’s orbit is not coplanar 
with the Earth’s orbit, the distances must be identical at a point where the object passes through the 
ecliptic, the plane in which the Earth orbits the Sun. Warning times are studied with the aid of coplanar 
orbits, and a condition required for a collision in the more general situation of noncoplanar orbits is 
examined. 
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Warning Time 

The amount of warning before a collision is the interval between the time the object is first discovered 
(or, to be more precise, the time the orbit becomes known accurately) and the time of collision. Although 
several orbital periods may elapse during this time, we concern ourselves here with the worst situation in 
which the comet or asteroid is discovered less than one orbital period before a collision. Unlike large 
asteroids with relatively short orbital periods, long-period comets (LPCs) do not present themselves for 
observation over multiple orbits, making it much more difficult to predict a collision decades in advance. 
Smaller asteroids may escape detection until less than one orbit remains before a collision. Warning time 
can be obtained through a straightforward application of time-of-flight equations. 

Figure 2 shows the Earth E in a circular orbit of radius r k = 1 astronomical unit (au) about the Sun S, 
and a comet (or asteroid) C in a coplanar elliptical orbit. The axes sq and .sq lie in the ecliptic, with sq in 
the direction of vernal equinox. It is assumed that discovery of C occurs after aphelion and before 
perihelion, at the point where the red line intersects the orbit of C, and E is assumed to be at the wrong 
place at the wrong time so that E collides with C before C reaches perihelion, at the point where the blue 
line intersects the orbits of both objects; the time of flight of C between these two points is to be 
determined. The motion of S and C is regarded as being governed by two-body orbit mechanics; that is, S 
and C are each treated as a particle, or a sphere whose mass is distributed uniformly, and the only forces 
exerted on 5 and C are those of mutual gravitation. The circular orbit of E about S also proceeds as two- 
body motion; however, the gravitational force exerted by E and C on each other is left out of account, as 
are all other perturbing forces acting on E, C, and S. 



Figure 2. Orbits of C and E. 
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The eccentric anomaly E at any point of an elliptical orbit obeys the relationship in equation (4.2-14) 
of reference 1 : 


cos E = 


1— r! a 


e 


( 1 ) 


where a and e are, respectively, the semimajor axis and the eccentricity of the orbit, and r is the distance 
between the centers of S and C at the point of interest. If E designates the eccentric anomaly at the point 
of discovery where the heliocentric distance is r, and E^ denotes the eccentric anomaly at the point of 
collision where the heliocentric distance is r k = 1 au, then the time of flight for an interval less than one 
orbital period is obtained via Kepler’s equation, as given by equation (4.2-9) of reference 1: 


t~t k = 



sin E)-{E k -e sinE^.)] 


( 2 ) 


where u is the gravitational parameter of the primary, in this case S. 

Now, the principal values of the inverse cosine function are 0 < cos _ l x < n; E is defined to be 0 at 
perihelion and n at aphelion. If E and E k are to have the correct values for the quadrants as illustrated in 
figure 2, the sign of the right hand member of equation (1) would have to be changed, as would the signs 
of E, Ek, sin E, and sin E k in equation (2); however, the absolute value of the time of flight t - t k would 
remain unaltered. 


Results 

Times of flight until collision for LPCs are calculated with the aid of equations (1) and (2). A family 
of orbits is constructed with perihelia r p equal to 0.1 and 1 au, and aphelia r a of 15, 20, 25, ..., 100, 
200, 300, ..., 1000, 2000, ..., 50 x 10 3 au. An aphelion of 15 au corresponds to an orbital period 
of about 20 years, whereas 50 x 10 3 au reaches the middle of the Oort cloud and corresponds to a period 
of 4 x 10 6 years. (LPCs are often regarded as those with periods greater than 200 years.) The distance r 
at which the orbit becomes known takes on the values 5, 6, and 7 au; corresponding time of flight, 
regarded as warning time, is presented in figure 3 as a function of r a . Warning time does not change 
appreciably for aphelia in the range 1000 < r a < 50 x 10 3 au, and a reduction in r p by a factor of 10 
reduces warning time by approximately the same amount as a reduction of 1 au in r. For the orbits 
studied, warning times range from 2 years in the case of a 20-year comet with r p = 1 au, detected at 
a distance of 7 au from the Sun, to 9.5 months in the case of a comet coming from the Oort cloud with 
r p = 0.1 au, detected at 5 au. 

Near-Earth asteroids (NEAs) are studied in a similar manner, with r p equal to 0.2 or 0.9 au, and 
r a = 1.4, 1.5, ..., 3.0 au; the associated orbital periods range between approximately 9 months, and 
2 years and 9 months. The detection distance takes on the values 1.1, 1.2, and 1.3 au, and warning time as 
a function of r a is presented in figure 4. Warning times range from approximately 90 days for an asteroid 
in a 0.9 x 1.4 au orbit, detected at a heliocentric distance of 1.3 au, to 7 days in the case of an asteroid 
with a 0.2 x 3.0 au orbit spotted when it is 1. 1 au from the Sun. 

Although it has been assumed thus far that the orbits of C and E are coplanar, a collision is possible 
also when the orbit plane of C is inclined to the ecliptic, so long as an ascending or descending node of C 
has a heliocentric distance of 1 au. Under these conditions, together with the two-body assumptions set 
forth in the discussion of warning time, the time of flight is unchanged. 
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Figure 3. Times of flight to collision, long-period comets. 



Figure 4. Times of flight to collision, near-Earth asteroids. 
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Designing a Collision 

Forthcoming batch filter analyses rest on the foregoing supposition, namely that a collision with Earth 
occurs when the object passes through the ecliptic plane and the object’s heliocentric distance is equal to 
r k , the radius at collision, or 1 au. By design, all of the orbits to be studied meet these conditions; this 
requires a certain relationship, to be developed presently, between the argument of periapsis on the one 
hand and, on the other, r p , r h and e. 

With the ecliptic chosen as the reference or fundamental plane, an object whose orbit plane is inclined 
to the ecliptic is, by definition, either at the ascending node or the descending node of the orbit when it 
passes through the ecliptic. Equations (1.5-4) and (1.5-7) of reference 1 allow us to write 

r= 'f 1 + e) (3) 

l + e cos v 

where r p is the radius of periapsis, e is the eccentricity of the orbit, and V is the true anomaly measured in 
the plane of orbit from the periapsis. At the ascending node v = -to, where co is known as the argument of 
periapsis, and at the descending node v = n - co. The requirement that r = r k at one of the nodes therefore 
can be expressed as 


r p (\ + e) 

1 ±e cos co 


(4) 


where the positive sign preceding cos co means the condition is imposed at the ascending node, and the 
negative sign is associated with the descending node. This relationship can be rearranged, 


±cos co = ■ 


1 + - - 


r k 


(5) 


Now, any member of equation (5) must of course remain between -1 and 1 (inclusive) if it is to be solved 
for co, 


-1 < ±cos co < 1 


( 6 ) 


however, one can simply work with the positive sign because it can be seen that the requirement is the 
same no matter which sign is used. Substitution from equation (5) into (6) gives 


( 


-1 < 


r k 


1 


1 


1 + - - -< 1 

V el e 


(7) 


or 


1 ZJL< r JL< x 
l + e r k 

and the right hand inequality yields the expected restriction on r p , namely 

r p ^ r k 


( 8 ) 

(9) 
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One can deal with the left hand inequality in (8) by appealing to equation (1.7-4) of reference 1 and 
substituting (r a - r p )/(r a + r p ) for e, which yields a rather sensible result involving the radius at apoapsis, 

r k * r a ( 10 ) 

Thus, as long as the inequalities (9) and (10) are satisfied, one is able to solve equation (5) for a value 
of argument of perihelion that meets the condition necessary for a collision to take place at the ascending 
or descending node, according to the choice of sign. (Changing the sign of the argument of the function 
cos -1 yields the supplementary angle: cos -1 ( x ) + cos -1 (-x) = n.) The analyst is presented with a second 
choice of sign because cos(co) = cos(-co); the result of cos _1 [cos(co)] provided by a calculating machine 
typically lies in the range of principal values of the inverse function, namely 0 < co < n, but a solution on 
the interval -n < co < 0 is also correct. Because co is always measured from the ascending node, co > 0 
always corresponds to a periapsis that is to the north side of the reference plane, and co < 0 always implies 
the periapsis is on the south side. Each of the four possible combinations of the two choices is associated 
with a pre- or postperihelion collision as indicated in table 1 . 

Table 1. Position of Collision 


Node of collision 

Perihelial hemisphere 

North (co > 0) 

South (co < 0) 

Ascending 

Preperihelion 

Postperihelion 

Descending 

Postperihelion 

Preperihelion 


Batch Filter for Long-Period Comets 

As mentioned in the introduction, the design of a system for performing observations can be guided by 
a quantitative study of the way in which orbit determination is influenced by the quality, number, and 
timing of the observations, as well as the number and location of the observatories. Because the orbital 
parameters of LPCs are distinct from those of asteroids, the investigation is divided accordingly with 
comets examined in the present section and asteroids in the following section. It would be highly 
desirable if a single system could serve effectively in observing the two classes of objects; however, it 
remains to be seen if this is in fact possible. 

In order to gauge the effects of the aforementioned factors on orbit determination, it is necessary to 
have in hand a way to judge the quality of orbital parameters obtained from observations. To this end the 
concept referred to as “erroneous predicted miss distance” is introduced. 

The study of LPCs begins with an introduction to the method of weighted least squares. The require- 
ment for collision expressed in equation (5) is used to construct a large number of orbital parameter sets, 
which then yield “true” values of measurements to which errors are introduced intentionally. The 
measurements thus corrupted account for limits in resolution that can be furnished by an actual telescope 
and are used to determine or estimate the orbit and obtain an associated erroneous predicted miss 
distance. The quality of preliminary orbit determination obtained from three observations is examined, 
followed by a study of the effects of multiple (more than three) observations and observations taken from 
two or more observatories. 
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Erroneous Predicted Miss Distance 


Orbital elements determined on the basis of observations generally differ from the true orbital 
elements, that is to say at a particular value of time t, the true position r and velocity v differ from the 
position r' and velocity v' resulting from the process of orbit determination. In practice the true orbital 
elements are unknown; however, they are to be specified in analysis that follows. The two sets of orbital 
elements can be compared in order to judge the quality of the result of orbit determination; however, it is 
more convenient to compare a single parameter if at all possible instead of six scalar values associated 
with position and velocity (or, for that matter, six classical orbital elements). Such a parameter e, 
hereafter referred to as an “erroneous predicted miss distance,” is now introduced. When e vanishes, the 
orbit determined from observations results in a collision at the specified position and at the designated 
time. (Other metrics can be found in the literature, but seem to require a comparison of two scalars.) 

An object C travels along a true orbit, shown in figure 5 with a solid red arc, designed to collide with 
the Earth E at time t^, at a point K where the object passes through the ecliptic plane at a heliocentric 
distance of r k = 1 au. The orbit of C is generally not identical to an orbit determined from measurements, 
shown with a dashed red arc, and associated with an object C' . An orbit solution is obtained for some 
epoch t, at which time the true position vector from S to C is r (t), and the position determined from 
observation is shown as r'(t). When the true orbit is propagated to the time of collision tk, r{tk) is of 
course the position vector from S to K. The orbit of C’ can also be propagated to A, yielding a position 
vector r '(tk); the magnitude £ of the difference r'(fy) - r(4) is defined to be the “erroneous predicted miss 
distance.” 



Figure 5. True orbit and orbit determined from observations. 
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There is one drawback to using e to judge the quality of an orbit solution: e may vanish even though 
r '(t) ^ r(t) and \'(t) ^ v(t). A solution to the Gauss problem, described in section 5.2 of reference 1, 
yields the velocities v(t 1 ) and v(t 2 ) at two specified positions r(t 1 ) and r (t 2 ), where the time of flight 
t 2 - /| between the two positions is also specified, as is the “direction of motion” (i.e., whether the orbit 
proceeds from r(tj) to r(/ 9 ) via an angular displacement less than or greater than n). Consequently, for 
any position solution r \t) and time of flight t /( - t, there exists some velocity v'(0 that results in C' 
passing through K at precisely t,, thus e = 0 even though the orbit of C' differs from that of C. If by 
unlucky coincidence an orbit determination algorithm happened to settle on just such a combination of 
r'(t) and v'(t), £ would vanish, giving a false impression that the orbit solution was identical to the true 
orbit. The effects of such coincidences will be diminished by introducing errors in the measurements in a 
random fashion, and averaging £ over a large number of measurement sets. 

For each orbit examined in the sequel, £ is obtained by the following steps: 

1. A set of classical orbital elements is constructed in accordance with equation (5) and used together 
with orbit propagation software to produce time histories of r and v for an interval of time over 
which optical measurements are to be taken. 

2. A time of flight to collision t, - t is determined according to equation (2), whereupon Lagrange 
coefficients F and G ( p. 179, ref. 2) are computed and used to obtain the position of collision r(L). 

3. Several values of time t t (i = 1, 2, 3, ...) are selected for obtaining observations. For each such 
value a position of an observatory is constructed and subtracted from r(/-) to yield the position 
vector from the observatory to the object of interest; this vector is in turn used to calculate a 
longitude and latitude of the object with respect to a set of inertially fixed heliocentric-ecliptic 
axes. 

4. For each observation, errors are intentionally introduced in the longitude and latitude to reflect the 
limits in angular resolution of a telescope; the resulting angles serve as measurements. 

5. A guess is made as to position and velocity at the time t of one of the observations. The object of 
the process of orbit determination is to adjust the guess and obtain r'(t) and v'(0 that are as 
consistent as possible with the measurements. 

6. The time of flight obtained in step 2 is used together with r \t) and \'(t) to obtain r'(t k ). 

7. The erroneous predicted miss distance £ is the magnitude of 8 = r \t k ) - r (t k ). 

Method of Weighted Least Squares 

Precision orbit determination can be performed with the linear method of weighted least squares 
formulated by Carl Friedrich Gauss in 1809, sketched out in references 3 and 4, and explained in detail in 
reference 5. We have used this method to obtain preliminary determination of orbits using three or four 
observations as described subsequently, as well as improved orbit determination from many observations 
and multiple observatories. 

The method of weighted least squares can be described briefly as follows: Let X\, be six orbital 
elements, for example, six scalars associated at a particular instant of time with position r(t) and velocity 
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v(0, to be determined from measurements of longitude and latitude 

obtained at the times t\, ...,t n of n observations. Using estimates Xj, ... ,X g of the orbital parameters, 

one may compute corresponding values of longitude (j^Xj, ...,X 6 ,qj, ...^X^ ...,X 6 ,t n J and latitude 

XfXj, . ..,X 6 ,?| ), ...,X(Xj, ...,X 6 ,t n ), and form 2 n residuals by subtracting the computed values from 


the measurements, 



A- \ 


^(x l5 .. 

■’*6 ’tj) 

(*1, 

..„x 6 ) = - 

Mtj-n) 

-X(x u 



) 


(;' = 1, ...,n) 
( j = n + 1, ...,2n) 


( 11 ) 


The purpose of the method of weighted least squares is to find the values Xj , ...,X^ that minimize the 
sum of the squares of the weighted residuals, 


2 n 

Q =y L w jy 2 j (i2) 

7=1 

where weights w\, ...,W 2 n are chosen to give more weight or importance to the measurements obtained 
with a better resolution; Wj is the assumed accuracy of the yth measurement. Now, for values of 

Xj, ...,X 6 in the neighborhood of Xj, . ,.,X 6 , the computed values can be expanded in a Taylor series, 
for example, 



which can be used to write approximate expressions for residuals 


y j{X\, ...,x 6 ) = $(tj)-<$>( 

X\, ■■■ ,X 6 ,t y) 





~yj 

(x u ... 

1 

In 

1 

50 

<>< 


% -’X 6 ,tj) 

(7=1, ■ ••,«) 

(14) 

If one defines arrays 







{y} = [yi{x h 

-,X 6 ) y 2 {x 1 , 

...,x 6 ) 

... y 2n (X 1 , ...,X 6 )] T 

(15) 

{y}± 

yi(^i. 

-A) y 2 (x i, 

...,x 6 ) 

y2n[X\i • ■ ■ > Xft ) 

1 T 

(16) 


A 


A \ 

/ - \1 T 




{x} = 

(Xi-Xj) (x 2 

-x 2 ) . 

■ (^6-^ 6 ) 


(17) 
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and matrices [H] and [W], whose elements are defined as 


Hi 


f (*■■ " 

■’^6 ,tj) 

(7=1, 

7=1, . 

..,6) 

g(*>- •• 


(j = n + l. 

. . . ,277 ; 7=1, 

..,6) 

W- 

yv Jl 

_ W i 
0 

O' = 7 ; 7=0 
0 * 7 ; 1=0 ••• 

,277.) 

,277) 



then Q can be expressed as 


2 n 


< 2=1 Wjyj 

j = 1 


= (l} T [W]{y} - ({y}- [H]{ x}) T [W]({y>- [H]{x}) 


( 18 ) 


(19) 


The value of {x} that minimizes Q is obtained by setting dQ/d{x) = {0}, a 6 x 1 array of zeros, which 
yields 

{x} = ([tf] T [W][//])”\//] T [W]{^ (21) 

as long as is nonsingular; this quantity is known as the normal or information matrix and is 

equal to the inverse of the covariance matrix describing the accuracy of {x} . When the normal matrix is 
singular, the orbit cannot be determined uniquely from the measurements and the orbit is said to be 
unobservable (ref. 5, sec. 4.12). A new estimate of the orbital parameters is formed by adding the 

T 

adjustment {x} to the old estimate (x} = [x 1 x 2 X3 x 4 x 5 xg] , and the process is repeated, leading to 

/v * /v * 

the values of x \ , . . . ,xg that minimize Q. 

As discussed in reference 5, each row of the 2 n x 6 mapping matrix [77] can be regarded as the product 
[77][<E»(7y,?)], where [ H ] is a 1 x 6 matrix of partial derivatives of a computed value ((]) or X) with 
respect to A | , . . . ,A 6 , and [0((y,t)] is the 6x6 state transition matrix for the time of the measurement tj and 

the time t at which the orbit is to be determined. Both matrices are evaluated with the values X\, ...,Xg. 
Although the rows of [ H ] are often computed with these products, it is also possible to form [ H ] 
numerically and this is one option available with the MATLAB® function LSQNONLIN. With this 

software, a search direction can be determined with the Levenberg-Marquardt method, a hybrid of 
the Gauss-Newton approach and the method of steepest descent (ref. 6, pp. 2-16 to 2-21). Parameter 
options that we set specifically for orbit determination include the minimum change in variables for 
finite differencing, DiffMinChange = 5 x 10 -4 , the termination tolerance on the function value, 

TolFun = 1 x 10 -10 , and Large Sc ale is set to of f to select the Levenberg-Marquardt method rather 
than a large-scale optimization algorithm. The elements of {x} are scaled so that the first three are in 
units of au and the last three are in units of km/s. The residuals are rectified so that -71 < yj < n and are 

then multiplied by 1 x 10 6 to avoid difficulties when the residuals are on the scale of the numerical 
precision of the computer, and where round off errors can cause a serious loss of accuracy. 
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Preliminary Orbit Determination 

The method of least squares has been used together with the seven-step procedure described in 
connection with the erroneous predicted miss distance, and MATLAB® software, to determine how the 
accuracy of preliminary orbit determination based on three observations is affected by the time interval 
between successive observations, by the distance r at which the object is detected, and by the resolution 
of the telescope. 

Proceeding with the first step, a set of classical orbital elements is specified as r a = 10 au, r p = 0.5 au, 
i = 16°, and Q = 36.5°, where i is the inclination of the object’s orbital plane to the ecliptic, and Q is the 
longitude of ascending node measured from the vernal equinox. The argument of perihelion co is 
computed according to equation (5) to be 93.017°. Observations are to be performed when the object’s 
heliocentric distance r is in the neighborhood of some specified value, so an initial value of true anomaly 
v(to) is obtained by solving equation (3) for v, chosen < 0 so that the object is traveling toward perihelion. 
With this set of orbital elements, time histories of r and v are recorded over an interval of 98 days. 

Continuing with step 3, the observatory is assumed to travel in the ecliptic plane in a circular orbit of 
radius 1 au; however, the observatory is not necessarily coincident with the Earth. At the beginning of 
the 98-day observation period, the observatory has an initial true longitude Lq = 180°, measured in the 
ecliptic from the direction of vernal equinox. The times of the three observations are assumed to have 
equal spacing, (* 3 - t 2 ) = {t 2 - t \ ), ranging from 5 to 49 days. (For example, measurements obtained at the 
beginning of the 1st, 50th, and 99th days have an interval of 49 days between successive observations, 
and the total span of the data arc is 98 days.) 

Figure 6 illustrates the comet’s orbit with a blue curve, and the red portion of the orbit indicates the 
interval over which observations are made, beginning at v(t Q ) = -159.7°. The green plane represents the 
ecliptic and the asterisk denotes the Sun. The three black points indicate the positions of the observatory 
at the beginning, middle, and end of the observation period, where the beginning is marked by the 
leftmost point. 

As described in step 4, each of the six measurements is regarded as the sum of a true angle and an 
error, which is now assumed to be no greater than the telescope resolution p. The MATLAB® function 
RAND is employed to produce pseudorandom errors uniformly distributed between the limits of the 
telescope’s resolution, ±p. In a Monte Carlo approach, an orbit is determined 100 times using a different 
random number seed in each trial, and an average erroneous predicted miss distance is reported. The 
assumed distribution of the error is not actually an important consideration when using the method of 
least squares since this method does not account for such statistics. A uniform distribution is simple to 
generate and just as valid as any other in this case. 

Figure 7 shows, on logarithmic scale of base 10, the average value of £ as a function of the interval 
between observations, (t 3 - t 2 ) or (t, - /, ), where the time of the second observation t 2 is taken to be the 
50th day of the 98-day data arc. The heliocentric distance r at C is varied between 5, 6, and 7 au by using 
initial values v(t Q ) = -155.3°, -159.7°, and -163.7°, respectively. Each point on a line represents the 
average value £ of £ in units of lunar distance from Earth (384400 km), from a set of 100 measurement 
errors formed with p = 0.2 seconds of arc (arcsec). Accuracy of orbit determination becomes better as 
r(t 2 ) decreases; however, the most dramatic improvements are obtained by increasing the observation 
interval to 49 days, at which point all three values of e are less than 2 lunar distances. Using the laws of 
logarithms, it is easily shown that a relationship of the form e = £o( ? 2 -? l) W > where e 0 is a constant, 
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Figure 6. Observation arc. 



Figure 7. Average £ for various detection distances. 
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Figure 8. Average £ for various angular resolutions. 


becomes log^o £ = m logio(t2 - h) + logio £o> which is the equation of a line with slope m. In each line 
of figure 7, £ decreases by two powers of 10 for every power of 10 increase in (t 2 - /,); therefore m = -2, 
meaning the accuracy of the preliminary orbit improves with the square of the observation interval. 

Figure 8 displays similar information, with r(t 2 ) held fixed at 6 au and p = 0.05, 0.10, and 0.20 arcsec. 
Telescope resolution is seen to have an effect on preliminary orbit determination accuracy, but the length 
of the observation interval is again the most significant factor. With an observation interval of 49 days, 
all values of £ are below 1 lunar distance. Unfortunately, longer observation intervals yield shorter 
warning times. 

In the same vein, an extensive analysis of preliminary determination of LPC orbits is carried out with 
the following additional details. 

In connection with the first step, 1008 sets of classical orbital elements are constructed with r a =10, 
40, 70, and 100 au, r p = 0.1, 0.4, 0.7, and 1.0 au, i = 10°, 30°, 50°, ..., 170°, and Q = 0°, 30°, 60°, ..., 
180°. Time of flight is nearly constant with respect to r a for r a > 100 au (see fig. 3); therefore, such orbits 
are not considered. For each set, the argument of perihelion co is computed according to equation (5). A 
98-day data arc begins with r(t Q ) = 6.5 au, and v(/ f) ) is chosen accordingly. 

As before, the observatory is assumed to travel in the ecliptic plane in a circular orbit of radius 1 au, 
and it is worth remembering that the observatory is not necessarily coincident with Earth. The effect of 
observatory position on orbit determination is studied by allowing the initial true longitude of the 
observatory to take on four values, each differing by 90°. Thus, a total of 4 x 1008 = 4032 cases are 
examined. The times of the three observations are equally spaced, with 33 days between subsequent 
observations. 
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Figure 9. Accuracy of comet orbits, as a function of i and r p . 

For step 4, a telescope resolution of p = 0. 1 arcsec is used to form 100 sets of measurement errors for 
each orbit to be examined, and the average value of erroneous predicted miss distance is recorded. 

The results are displayed in figure 9 as functions of i and r p \ each data point represents £ for one case, 
and a family of 112 cases (in which r a , Q, and vary) form one stack of points. The best orbit 
determination is obtained when the orbital planes of the observatory and comet are perpendicular, 
whereas nearly coplanar orbits yield the poorest determinations. This relationship to inclination stems 
from the fact that, with three observations, an orbit is unobservable when it is coplanar with the orbit of 
the observatory. One may form analytic expressions for the elements of the matrix [if] discussed in 
connection with the method of weighted least squares which are partial derivatives of longitude and 
latitude. In addition, one may determine which elements of the state transition matrix [<h] are zero when 
two-body mechanics are assumed; a reference orbit coplanar with that of the observatory is used and [<1>] 
is written according to equations (9.84)-(9.87) of reference 2. It can then be shown that for three 
observations, the 6x6 matrix [H] is singular, but for four observations the 6x6 matrix [//] 1 [H] is not 
singular. A study of figure 9 shows that with r p = 1 au, retrograde orbits are harder to determine 
accurately than prograde orbits. Orbit determination accuracy with r p = 0.1 au is noticeably poorer than 
with r p = 0.4 or 0.7 au, and even poorer with r p = 1 au. In figure 10 the results are shown as functions of 
r a and r p ; each stack contains 252 values of £ associated with varied i, Q, and L {) . Orbit determination is 
the most difficult when r a has the smallest value, 10 au. 

Results are not presented as functions of Q, co, and v(/q) because there is limited value in doing so. No 
remarkable relationship appears to exist between £ and Q. The angle v(/ () ) simply determines the point of 
the object’s orbit at which observations begin. Finally, argument of perihelion co is dependent upon r a 
and r p through equation (5). 
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Figure 10. Accuracy of comet orbits, as a function of r a and r p . 


The 4032 average values of £ are sorted from largest to smallest, and the worst and best cases 
are identified. The effect of the length of observation interval is then studied for these two cases, 
with times of the three observations again given equal spacing, ranging from 5 to 49 days between subse- 
quent observations. The worst case, numbered 246, is associated with the orbital parameters r a = 10 au, 
r p = 1 au, i = 170°, Q = 0°, co = 0°, and v(t Q ) = -151.68°, and Lq = 180° for the observatory. Curves for 
the cases with the same orbital parameters and other values of Z, Q are also shown in figure 11. 
Observatory location affects the accuracy of orbit determination, but not to the same extent as length of 
the data arc. In all four cases, £ is less than 4 lunar distances with an observation interval of 49 days. 

The best case is numbered 3946, with orbital parameters of r a = 100 au, r p = 0.7 au, i = 110°, 
Q = 120°, co = 66.69°, and v(7 0 ) = -142.86°, and L Q = 90° for the observatory. Curves for the cases with 
the same orbital parameters and other values of Lq are shown in figure 12; long observation intervals once 
again result in the best preliminary orbit determination. With an observation interval of 49 days, the 
value of £ associated with each observatory is less than 0.1 lunar distance, or about 6 times the Earth’s 
radius. A slope of m = -2 is evident in figures 11 and 12, confirming the inverse square relationship 
between £ and observation interval established in figures 7 and 8. 

The foregoing results involving preliminary orbit determination corroborate statements made in 
references 3 and 4, pointing out that the length of the data arc is the single most important factor in 
determining the accuracy of the orbit solution. The number and precision of the observations, the object’s 
proximity to the observatory when measurements are obtained, and even the use of radar measurements 
are all secondary to the length of the data arc. 
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Figure 11. Worst orbit, various observatory longitudes. 
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Figure 12. Best orbit, various observatory longitudes. 
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Multiple Observations 

As opposed to the classical methods of Laplace, Gauss, and Olbers, the method of least squares 
permits the use of a great number of measurements, which can be obtained from one or more 
observatories. We now proceed to study the effects of multiple observations and observatories on orbit 
determination of LPCs. A number of observations, ranging from 3 to 99, are taken in equally timed 
increments over a period of 98 days. We have chosen to obtain an orbit solution for the time of the first 
observation. Differential correction begins with an estimate {L} corresponding to an arbitrary day, 
typically a day after the end of the 98-day observation period. 

The forthcoming results, shown in figures 13 through 18, are reported first for a single observatory, 
and then for multiple observatories. 

Single observatory. The orbits identified in the discussion of preliminary orbit determination as result- 
ing in the largest and smallest average value £ of e, cases 246 and 3946, respectively, are revisited. 
Figures 13 and 15 show results calculated for various numbers of observations obtained from a single 
observatory. In each case, a solid curve shows the average erroneous predicted miss distance £ from 
100 sets of measurements, each of which is produced with a telescope resolution p of 0.1 arcsec. 

Results for the worst orbit are contained in figure 13. With 3 observations spaced 49 days apart, e is 
nearly 3.8 lunar distances, whereas 99 observations spaced 1 day apart reduce £ to less than 0.22 of a 
lunar distance, or about 13 Earth radii. In the case of the best orbit, the quality of the result is evidenced 
by the use of Earth radius (6378 km) as the unit for measuring average erroneous predicted miss distance. 
The data in figure 15 show that with 3 observations e is approximately 2.5 Earth radii; 99 observations 
improve the measure to less than 0.5 Earth radius. In the worst case there is a marked improvement from 
increasing the number of observations from 3 to 5, followed by a more gradual slope of approximately 
m = -1/2 in going from 10 to 99 observations. In the best case, the slope is also approximately -1/2. As 
discussed previously, this is an indication that the accuracy of the orbit as measured by £ improves as the 
square root of the number of observations. 

As one might expect, better telescope resolution (indicated by a smaller value of p) improves orbit 
determination accuracy. In order to quantify the improvement available with very high resolution 
measurements (obtained, for example, with interferometry) an analysis is made with 1 1 observations 
taken over a 98-day data arc involving a mixture of resolutions; p = 0.01 arcsec for initial observations, 

and p = 0.0001 arcsec for some number of final observations. As described earlier, weights are assigned 
using the value of p with which the corresponding measurements are made, thus giving more weight or 

9 

importance to the measurements obtained with a better resolution (wj = I / p / )- In figures 17 and 18, 

improvements in £ obtained for the worst orbit (Case 246) and the best orbit (Case 3946), respectively, 
are shown. The number of final observations made with the better resolution is indicated on the abscissa. 
A comparison of figures 13 and 17 reveals £ is made better by a factor of 10 when the resolution of 11 
observations is improved from 0.1 to 0.01 arcsec. The same conclusion is reached by comparing fig- 
ures 15 and 18. In each case, a further improvement by a factor of 10 is obtained when the last 
4 observations are made with p = 0.0001 arcsec. 
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Figure 13. Worst orbit, various observatory configurations. 



Figure 14. Worst orbit, various observatory configurations (detailed view). 
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Figure 15. Best orbit, various observatory configurations. 



Figure 16. Best orbit, various observatory configurations (detailed view). 
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Figure 17. Worst orbit, observations with mixed resolutions. 
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Figure 18. Best orbit, observations with mixed resolutions. 
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Multiple observatories. As previously mentioned, the method of least squares permits orbit determina- 
tion to take advantage of measurements not only from additional observations, but also from additional 
observatories. The benefits of configurations of 2, 3, and 4 observatories are studied with the aid of the 
worst and best LPC orbits. 

The positions of two observatories are constructed such that they have heliocentric circular orbits of 
radius 1 au in the ecliptic plane, and their true longitudes L Q are always 180° apart. Three observatories 
are placed in similar orbits, with their true longitudes phased by 120°. In the case of four observatories, 
the first two have orbits identical to the configuration of two observatories just described, and the 
remaining two are in similar coplanar orbits perpendicular to the ecliptic with the ascending node Q 
taking on values of 0°, 45°, and 90° in order to determine what effect, if any, Q has on orbit 
determination. 

In connection with the worst comet orbit, figures 13 and 14 display e as a function of the number of 
observations for all observatory configurations described heretofore. It is evident that two or more 
observatories offer an improvement in e of nearly a factor of 10 over that from a single observatory. 
With 99 observations each, two observatories yield 8 less than 0.016 lunar distance, slightly less than 
1 Earth radius. Three observatories are only marginally better than two observatories. Four observatories 
are not substantially better than three, and 8 is relatively insensitive to the value of Q for the members of 
the four-observatory configuration that have orbits normal to the ecliptic. Similarly, figures 15 and 16 
show that in the best case substantial improvement in 8 is obtained by employing two observatories 
phased by 180° instead of a single observatory, but the addition of a third or fourth observatory does not 
appear to be cost effective for this orbit and observatory heliocentric distances of 1 au. 

The advantage of two observatories over one may stem from a doubling in the number of observa- 
tions, or from parallax; in order to gauge the contribution of each, one may form a configuration of two 
observatories that are coincident with one another, thus eliminating parallax. Results obtained from this 
configuration are represented with a dash-dot curve in figures 13 and 15. The first of these shows that in 
the worst case, the improvement obtained from a pair of observatories stems more from parallax than 
from doubling the number of observations. Figure 15 shows that in the best case, the benefit of two 
observatories has as much to do with parallax as it does with the availability of twice as many 
observations. 

Of the multiobservatory configurations examined here, two observatories in circular heliocentric orbits 
of radius 1 au in the ecliptic plane, phased 180° apart, provide the best balance of cost and accuracy in 
determining the orbits of LPCs. Thirty to forty observations equally spaced over a 98-day interval appear 
to give results nearly as good as 99 observations taken 1 day apart. 

Sequential Filter for Long-Period Comets 

The method of weighted least squares, described previously and used to produce all of the results 
presented earlier, is referred to as a batch processor because it employs observations collected together in 
a single set; the larger the batch, the more arithmetical operations are required to carry out the matrix 
multiplications indicated in equation (21). Moreover, when new observations are obtained, they must be 
added to the previous batch to form an even larger set, and the entire orbit determination procedure must 
be repeated. Another shortcoming of the method is that there are no provisions for treating the state and 
observation deviations as random processes and using any associated statistical information to improve 
the orbit solution. 
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Remedies for these problems exist in the Kalman filter, a recursive or sequential algorithm in which 
measurements are used continuously, as they become available, to refine the solution for the orbit. The 
recursive property of the procedure allows the previous solution to be used as the starting point for the 
solution to follow. The computational expense of the sequential filter can be less than that of the batch 
filter because the matrix requiring inversion in the former case can be smaller — when only one 
measurement is processed at a time, the matrix is reduced to a scalar. 

We describe briefly two variations of the Kalman filter, the conventional and the extended algorithms. 
A method for obtaining the a priori information needed to start the filter is also discussed. Warning time 
offered by a sequential filter with high resolution optical measurements from a single observatory is 
compared with that derived from measurements of lower resolution from two observatories. The advan- 
tages of measurements of range and range-rate are illustrated. Effects of comet outgassing are con- 
sidered, followed by a discussion of the construction and use of the body-targeting plane or B-plane in 
orbit determination analysis. Finally, the calculation of probability and probability of collision are taken 
up. 

The Kalman Filter 

The conventional Kalman filter is developed from differential equations for the state, and expressions 
that relate the measurements to the state, both of which have been linearized about a reference state time 
history. To reduce the effects of the errors introduced in the process of linearization, the reference state 
can be replaced after each observation with the improved estimate of the state; this procedure is referred 
to as the extended Kalman filter. The two varieties of the filter are discussed in turn in what follows. 

Conventional Kalman filter. The conventional Kalman filter is described in detail in section 4.7.1 of 
reference 5 in terms similar to those used earlier in connection with the batch algorithm. 

The n x 1 column matrix {X} contains n {n = 6) orbital elements, dot products of the position r(7) and 
velocity \{t) with three right-handed, mutually perpendicular unit vectors. The state deviation matrix \x ) 
is defined as the difference between the true elements {X} and the elements { X *} associated with a 
reference orbit, 




( 22 ) 


The m x 1 observation deviation, or residual matrix {y}, is defined as the difference between m 
measurements {T} and the values {7*} calculated from the reference orbit, 

{y}={Y}-{Y*} (23) 


where the measurements {7} are regarded as the sum of nonlinear functions {G({X}, /)} of the state and 
time, and errors {ej, 


{7} = {G({X},t)} + M (24) 

Equations for the filter time update are given as equations (4.7.18) and (4.9.50) of reference 5: 

{Xj} = \<&( tj ( 25 ) 
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(26) 


+[Qi- 1 ] 




lj-\\ 


where {xj_ 1 } and {x j } are, respectively, the estimated values of {x} following a measurement update at 
an instant of time tj- \, and prior to a measurement update at tj. The n x n state transition matrix 
1 d>( t : , t :_ | )] that relates the state deviations at times tj- \ and tj is discussed presently. The n x n estimate 
error variance-covariance matrices (or simply covariance matrices) of the state deviations {*j- \ } and 
{Xj} are denoted by [Pj-\] and [Pj], respectively. The n x n matrix Qj~\ is the covariance at tj- \ of a 
zero-mean, white sequence known as process noise or state noise. 

The m x 1 observation deviation {yj}, m x n observation-state mapping matrix [H j], and the n x m 
Kalman gain [K j\ are assembled in preparation for the measurement update: 



(27) 


(28) 


[^] = [P,][H;] T ([H;][P j ][H ;] T +[^ 1 


(29) 


where { Yj} contains the actual measurements at tj, and where the m x m matrix [7?y] is the covariance at tj 
of {e}, assumed to have a normal distribution and a mean of zero. Expressions for the elements of [// ,■] 

associated with measurements of longitude and latitude are developed presently, and measurements of 
range and range-rate are taken up subsequently. 

Equations for the filter measurement update are set forth in equations (4.7.16) and (4.7.12) of 
reference 5, 


{x j } = {x j } + [K j ]({y j }-[H j ]{x j }) 

(30) 

[P j ] = ([I]-[K j ][H j ])[Pj] 

(31) 


where {xj } contains the estimated values of the state deviation following a measurement update at time 
tj, and [Pj\ is the corresponding covariance matrix. After any update performed with measurements 
obtained at tj, the best estimate (A' } of the orbital elements can be formed by adding {xj } to the 

elements { Xj } associated with the reference orbit at that time, { X'j } = { Xj} + { Xj } . 


In the forthcoming discussion of Gauss’s method we take up the matter of obtaining the a priori 
information needed to start the filtering process, namely the initial values (Aq } of the elements of the 
reference orbit, the initial state deviation estimate {v 0 }, and a consistent initial covariance matrix [Pq], 
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The state transition matrix \0(tj,tj_ i)] required in equations (25) and (26) for the time update is 
formed on the assumption that motion of the Sun and the object is governed by two-body gravitational 
mechanics. The position vector r*_j at time tj- \ from the mass center of the Sun to the mass center of 

the object traveling on the reference orbit, and the inertial time derivative v (_ ] of that position vector, are 
related to the relative position rj and velocity v* at another time tj by equations (9.68) of reference 2, 

r] = Fr]_ ] +Gv*_ ] (32) 


* 



j-f * 

= F t r i- 1 


+G ^U 


(33) 


where F, G, F t , and G t are known as Lagrangian coefficients, which can be calculated with Battin’s 
universal variables as indicated in equations (9.69) of reference 2. The state transition matrix is 
partitioned into four 3x3 parts, 


Wfj'tj-i)} 


[R(tj,tj_ i)] [R(tj,tj_ i)] 
o] [v(tj,tj_ i)] 


(34) 


where the superscript * is omitted because it is used by Battin to indicate the adjoint of a matrix, but it is 
to be understood that the partitions are evaluated with r*_i, v*_j, r*, and v* of the reference orbit 
according to equations (9.84)— (9.87) of reference 2. The universal variables U 0 , Lj , U 2 , and U 2 required 
to evaluate the Lagrangian coefficients are obtained straightforwardly from the relationships given in 
problem 4-21 of reference 2. The coefficient C, required to form [V(tj,t ,_i)] according to equa- 
tion (9.87), is expressed in equation (9.74) in terms of U 2 , t/ 4 , and U 5 . The second of these, ( 7 4 , is easily 
obtained from C/j, U 2 , and U 2 with the aid of equation (4.108). The universal variable U 5 can be 
expressed as in problem 4-30 in terms of t/j and U 2 , except when % = IT) = 0, in which case U 0 = 1, and 
U 2 = U 2 = C/ 4 = 0; the variable u also vanishes in view of equation (4.100), which leads to q = 0 by way 
of equation (4.104), and to 6/5(2%) = 0 from equation (4.1 12). Finally, the first of equations (4.1 13) yields 

6/5 (x) = W) = 0. 

In connection with measurements of longitude § and latitude X, the observation-state mapping matrix 
[H A defined in equation (28) is derived as follows. The position vector from the Sun S to the object C 
can be written as 


r- r‘ SC - XjSj + X 2 s 2 + X 3 S 3 (35) 

where sj, S2, and S3 are a set of right-handed, mutually orthogonal unit vectors fixed in an inertial or 
Newtonian reference frame N. Unit vectors s | and §2 lie in the ecliptic plane, with s | in the direction of 
vernal equinox, and the direction of S3 is north of the ecliptic. Similarly, the position vector from S to an 
observatory O can be expressed as 


r‘ so - OjSj + 0 2 s 2 + O3S3 


(36) 
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and the position vector d from O to C is thus 


d= r oc = r 5 C -r 50 =(X 1 -0 1 )s 1 +(X 2 -02)s 2 + (X 3 -0 3 )s3 


(37) 


A unit vector d that has the same direction as d can be brought into a general orientation in N by first 
giving it the same direction as s 3 and then subjecting it to a body-two, 3-2-3 rotation sequence with 

angles of <]) (longitude), -A (negative of latitude), and 0 (zero). Consequently, d can be expressed as 


d = cos (|) cos As j + sine]) cos As 2 + sinAs 3 (38) 

and relationships for longitude and latitude follow, 


d So d So 

tancj) = — — — = — -X 

(39) 

d Sj d s 3 


sinA = ds 3 = 83 

3 d 

(40) 


1 / 2 

where d is the magnitude of d, d = (d d) 


Each observation of an object with an optical telescope yields measurements of (|) and A; the 
corresponding matrix | H ] ■ is then 
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(41) 


The partial derivatives of <j) are obtained by noting that 
equation (39), 


— tan 1 u = — — ; therefore, in view of 

dx 1 + u 2 dx 


9 (() _ (d-Sj ) 2 3 (d-s 2 ) 

dX{ (d-s 1 ) 2 + (d-s 2 ) 2 dX t [(d-sj) 


(d-sQ 2 

(d-s 1 ) 2 + (d-s 2 ) 2 


_L dd_ . 

(d-s^axT* 2 


(d-s 2 ) 3d . 

(d-s ^ 2 dX, 


(d-s 1 ) 2 + (d-s 2 ) 2 


~ dd „ » 3d . 

(d-S 1 ) — .^-(d-s,) — -S 1 


0 = 1 , ..., 6 ) 


(42) 
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With the aid of equation (37), it can be seen that 9d/9X ( - = s, for i = 1, 2, 3, and 9d /dXj = 0 for i = 4, 5, 6; 
thus, 


jj> _ ~(d-s 2 ) 

a ^i (d-s 1 ) 2 + (d-s 2 ) 2 

jj> = (d-sQ 

dX 2 (d-s 1 ) 2 +(d-s 2 ) 2 

■a.. 

9(|) _ _ 9(f) 

c)X 4 ~ dX 5 ~ dX 6 ~ 


(43) 


(44) 


(45) 

(46) 


The partial derivatives of A are obtained by noting that — sin 1 u = 


9x 


7 


1 — u 


du 

2 9x 


; therefore, in view of 


equation (40), 
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2d 3 dXj 


1 

9d „ (d-s 3 ) J 9d 

A /(d-S 1 ) 2 + (d-S 2 ) 2 

[ax, 53 a 2 ®, J 


( i = 1 , ..., 6 ) 


(47) 


from which we obtain 

9A _ -(d-s 1 )(d-s 3 ) 

9 *i 9 2 A /(d-S 1 ) 2 +(d-§ 2 ) 2 

9 A, = -(d-s 2 )(d-s 3 ) 

9X 2 t/ 2 -J(d-s 1 ) 2 +(d-s 2 ) 2 

9A _ ^(d-S}) - + (d-s 2 )~ 
dx^~ 

9A _ 9A _ 9A _ 

c)X 4 ~ dX 5 ~ dX 6 ~ 


(48) 


(49) 


(50) 

(51) 
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Extended Kalman filter. The extended form of the Kalman filter is presented in section 4.7.2 of refer- 
ence 5; as mentioned previously, it serves to reduce the effects of errors resulting from linearization upon 
which the conventional filter rests. The error reduction leads to faster convergence of the extended filter 
in comparison with the conventional form and is accomplished by replacing the reference orbit with the 
current, best estimate of the true orbit after each observation. As a consequence, the estimate of 

the state deviation following a measurement update vanishes, and it becomes unnecessary to perform the 
time update of the state deviation indicated in equation (25). The time update of the covariance, equa- 
tion (26), and preparation for the measurement update, equations (27)-(29), proceed as before. It can be 
seen that, with {xj } = {0}, the measurement update for the state deviation, equation (30), gives way to 

{x j } = [K j ]{y j } (52) 

The measurement update of the covariance, equation (31), remains the same, and the aforementioned 
replacement can be expressed as 


{ X* } ^ { X* } + [Kj]{ yj } 

(53) 

{*;}<-{ 0} 

(54) 


A Priori Information Obtained by the Method of Gauss 

The venerable method of Gauss for preliminary determination of orbits from optical measurements has 
been in use for two centuries; the occasion for its creation was the discovery (and subsequent disappear- 
ance behind the Sun) in 1801 of the first minor world, Ceres. The asteroid was recovered through the 
application of the Gaussian method with observations that spanned only 1 month, and Carl Friedrich 
Gauss became recognized immediately as the premier mathematician in all of Europe. 

The method requires only three observations consisting of two angular measurements each, which is in 
most cases the minimum number of observations needed. (As discussed previously, four observations are 
required when the orbits of the object and observatory are coplanar.) Gauss’s method is used to obtain 
the a priori information needed to start the filtering sequence, as described in what follows. 

A variant of Gauss’s method presented (but not labeled as such) in section 5.8 of reference 1 is 
brought to bear, with the Sun substituted for Earth in the role of primary body, and ecliptic longitude and 
latitude used in place of right ascension and declination. Battin’s expressions for the Lagrange 
coefficients F and G in terms of universal functions (eqs. (4.84), ref. 2) are used in place of the infinite 
series suggested in reference 1. The heart of the procedure is an iterative solution of a system of six linear 
algebraic equations, yielding estimates of the position r'(f 2 ) and velocity v'(f 2 ) at the time f 2 of the 
second of three observations. Such a solution is obtained with these steps: 

1. Three pairs of measurements of longitude (j), and latitude Xj for times /,■ (/ = 1, 2, 3) are produced 
according to steps 1, 3, and 4 of the procedure given in the discussion of erroneous predicted miss 
distance. The intentional errors to be added ({^},- in eq. (24)) are produced with the MATLAB® 

function RANDN, yielding a normal distribution with zero mean (in accordance with the Kalman 
filter assumption) and a standard deviation of 1. Errors are then multiplied by the telescope 
resolution p, but subsequently limited in magnitude; errors less than -2p are replaced with -2p, 

and errors greater than 2p are replaced with 2p. 
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2. The measurements constructed in step 1 are used to form three unit vectors d, with the aid of 
equation (38). 

o/t ^ 

3. An initial guess for the vector r'(f 2 ) = r (t 2 ) + d d 2 is formed with the kn own observatory 

SO 

position r (t 2 ) and a guess for the value of d, such as 5 au. 

4. An initial guess of the velocity v'(t 2 ) is produced by assuming that the object is in a circular orbit 
and therefore has a magnitude of ^/p/r 2 , where r 2 is the magnitude of r'(f 2 ). The direction of 
v'(f 2 ) is given by a unit vector normal to r '(t 2 ), having the same direction as (d| x d3)x r'(f 2 ), 
where the cross product dj x d 3 is approximately normal to the orbital plane. 

5. Lagrange coefficients F\ and G\ that relate r \t{) to r'(f 2 ) and v'(f 2 ) are constructed from these 
two vectors and the time difference t\ - t 2 . Coefficients F 2 and G3 corresponding to t 3 are 
obtained similarly. 

6. An iterative solution of equations (5.8-10) in reference 1 is terminated when the magnitudes of all 
six components of r \t 2 ) and v'(f 2 ) change by less than 1 x 10 -14 percent from one pass to the 
next, or after 50 passes, whichever occurs first. It is asserted in reference 1 that this process 
converges quickly if the time intervals t 3 - t 2 and t 2 - /| are “not too large.” 

7. Initial elements { Xq } for a reference orbit are simply dot products of r'(f 2 ) and v'(h) with S3, 
s 2 , and S3. Initial values of a state deviation {x 0 } are given by dot products of the same unit 
vectors with the vector differences r'(t 2 )-r(t 2 ) and x'(t 2 )-y(t 2 ), where r(l 2 ) and v(f 2 ) are the 
true position and velocity at t 2 . 

This procedure is used 100 times with 100 sets of observations to obtain 100 different matrices {Xq} 
and {vq}- Convergence was achieved in all cases. An initial covariance matrix [F Q ] describing the 
uncertainties in {v 0 } is calculated from the 100 state deviation matrices using the MATLAB® function 
COV. This a priori information is used to start a Kalman filter sequence. When Monte Carlo simulations 
are performed, each of the 100 pairs of associated matrices {Xq} and {x 0 } are used in a separate 
simulation, but the same matrix [P 0 ] is used in each simulation. 

Single Observatory, and Two Observatories 

The results presented in connection with multiple observations indicate that it may be possible to make 
reliable forecasts with two observatories whose angular resolutions are on the order of 0.1 arcsec, or with 
a single observatory whose resolution is better by 1 to 3 orders of magnitude. An extended Kalman filter 
is used to compare the amount of warning time offered by the two alternatives, where warning time is 
defined (somewhat more rigorously than before) to be the interval between the time a filter yields an 
erroneous predicted miss distance less than 1 Earth radius, and the designed time of collision. The task is 
facilitated by introducing a new term, “watch time,” defined to be the interval between the time £ is 
calculated to be less than or equal to 1 lunar distance, and the designed time of collision. A watch is 
issued prior to a warning for hurricanes and other reasonably predictable natural calamities, and so it 
would be in the case of a celestial collision. 
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Single observatory. We assume the single observatory achieves a resolution of 0.0001 arcsec by means 
of interferometry; two telescopes, each having a resolution of 0.01 arcsec, are placed some distance apart 
in a circular heliocentric orbit in the ecliptic plane with a radius of 1 au. Observations are obtained at 
7-day intervals from only one of the telescopes until e < 1 lunar distance, and the watch time is calculated. 
Subsequently, both telescopes are used together to produce measurements at 1-day intervals with 
p = 0.0001 arcsec, and the filter continues to process those measurements until e<l Earth radius, at 
which point a warning time is computed. 

Two observatories. The configuration of two observatories is the same as that previously described: 
two telescopes are placed in a circular heliocentric orbit in the ecliptic plane, with a radius of 1 au, and 
phased 180° apart. For the purpose of making the comparison, one of the telescopes is assumed to be 
coincident with one of the instalments in the single-observatory configuration. Each telescope has a 
resolution of 0.01 arcsec. For any one observation, this configuration produces twice as many 
measurements as the single observatory. Furthermore, two observatories provide parallax that cannot be 
obtained with the single observatory. Measurements are taken at 7-day intervals until e < 1 lunar 
distance, at which time the frequency increases to once per day until e < 1 Earth radius. 

All measurements contain Gaussian noise and are constructed as described in step 1 of the procedure 
for obtaining a priori information. 

Results. Watch and warning times have been computed for both configurations, using the 1008 
hypothetical LPC orbits given previously, and four values of the observatory initial true longitude L {) . 
The results obtained for L 0 = 180° are presented in figure 19. Computational expense prevents us from 
performing Monte Carlo analysis and presenting averages; each data point represents a single set of 
random measurement errors for that orbit. 

In almost all cases two observatories, enjoying the advantages of parallax and twice as many 
measurements, provide a greater watch time than a single observatory taking measurements at the same 
resolution. The advantage passes to the single observatory after it begins to obtain interferometric 
measurements; the improvement in resolution by a factor of 100 quickly reduces e to 1 Earth radius or 
less, giving a greater warning time in general than two observatories. 

Sixteen groups can be identified in the watch and warning times given in figure 19. The first, second, 
third, and fourth sets of four groups correspond to r a = 10, 40, 70, and 100 au, respectively. The four 
groups within each set are associated with r p = 0.1, 0.4, 0.7, and 1.0 au. The results for other values of L 0 
are similar to those shown in figure 1 9. 

When a conventional filter is used in place of the extended algorithm, it ceases to converge on the 
correct orbit and yields a warning time of zero in approximately 100 of the 1008 orbits, rendering it 
unsuitable for obtaining watch and warning times. An example of this behavior is exhibited in figure 20. 
The lack of an update to the reference trajectory in a conventional filter is the source of this defect. 

Range and Range-Rate Measurements 

The accuracy of orbit determination can be improved greatly with measurements of range and range- 
rate obtained from radar or lidar instruments. Range is simply the distance d between the observatory and 
the object, or the magnitude of the vector d given by equation (37). Range-rate, d, is the time derivative 
of d. Measurements are performed by analyzing the returning echo from an electromagnetic pulse that 
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Figure 19. Watch and warning times. 



Figure 20. Lack of convergence with a conventional filter. 
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has been directed toward the object. The distance d is related to the wave’s round trip time, and d is 
related to the Doppler shift in the frequency. 

As pointed out in the discussion of the Kalman filter, the observation-state mapping matrix [//.] is 
dimensioned m x n , the measurement covariance matrix [Rf] is in x m, and the observation deviation 
matrix {}>■} is m X 1, where m is the number of measurements and n is the size of the state, in our case 

n = 6. So far, only angular measurements of longitude and latitude have been considered, corresponding 
to m = 2. The addition of measurements of range and range-rate increases the value of m to 3 and 4, 
respectively. The covariance matrix [7?y] is diagonal; with m = 4, the first two nonzero elements of [Rj\ 
contain the square of the resolution of the optical telescope, the third is the square of the resolution of 
range measurements, and the fourth is the square of the resolution of range-rate. When observing an 
asteroid at 20 times the lunar distance, current terrestrial radar systems can achieve resolutions of 10 m 
and 0. 1 mm/s for range and range-rate, respectively. 


Expressions for the elements of [H j] associated with measurements of longitude and latitude are 
developed in equations (43)— (46) and (48)— (5 1). The third and fourth rows of [H j], associated respec- 
tively with measurements of range and range-rate, are derived as follows: 


dd 

dXj 


s | (dd) " 2 
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As mentioned earlier, it is evident from equation (37) that dd /d A, = s, for i = 1, 2, 3, and dd fdX l = 0 for 
i = 4, 5, 6; therefore, 


dd _ d s ( - 
dX f ~ d 



(i = 1, 2, 3) 


(i = 4, 5, 6) 


(56) 
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Range-rate can be expressed as 


d = — [(d- d) 1/2 ] = -(d- dr 1/2 2d- 

dt 2 


dt d dt 


(58) 


where N dd/dt indicates differentiation of d with respect to time in N, a Newtonian reference frame in 
which unit vectors s, (i = 1, 2, 3) are fixed. Strictly speaking, specification of a reference frame is not 
required because the left hand member of the equation is the derivative of a scalar, but the choice of N 
leads right away from equation (37) to 


N d ■ ■ ■ 

— — d= (X\ -Oj)Si + (^2 ~ 6>2)S2 _ 6?3)S3 

dt 


(59) 
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Now, 
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The last three elements of {X} are simply 
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therefore, in view of equation (59), — — ( N d d /dt) = 0 for i = 1, 2, 3, and — — ( N d d /dt) = s ; _3 for i = 4, 5, 6. 

oXj dXj 

Together with the values for 3d/3X ; - used previously, this leads to 
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The benefits of measurements of range and range-rate are illustrated in figure 21, with curves showing 
a reduction in average erroneous predicted miss distance 8 in units of Earth radius (R/;) as a function of 
the number of observations, taken once a day. The average is taken from 100 extended Kalman filter 
simulations, each with a different set of normally distributed measurement errors, and a priori information 
as discussed earlier. The analysis involves a hypothetical LPC designed to collide with Earth, with orbital 
elements obtained as described for designing a collision: r a = 100 au, r p = 0.7 au, i = 50°, 

Q = 60°, co = 66.68597°, and v 0 = -137.06483° (corresponding to r = 5 au). The orbital period of this 
comet is 1839.4 years. 

All measurements are obtained at a single observatory. The least accurate determination of the orbit, 
shown with the solid line, is obtained from measurements that are strictly angular with a resolution of 
0.01 arcsec. The slope of the line is approximately -2.3; therefore 8 varies inversely with the number 
of observations to the power 2.3. Greater accuracy is obtained when the angular measurements are 
supplemented with range measurements having a resolution of 10 3 km, as shown with the dashed curve, 
and the inclusion of range-rate measurements with a resolution of 1 m/s leads to the most accurate 
orbit determination, as displayed with the dash-dot curve. The usefulness of range and range-rate 
measurements is most pronounced over a short data arc of 15 to 30 days; they do not appear to be 
necessary over an arc of 80 or 90 days. 

Modeling Comet Outgassing 

A visible comet in the night sky typically appears as a fuzzy object; light is reflected diffusely from a 
cloud, or coma, consisting of dust evaporated by solar heating from the comet’s core, or nucleus. Solar 
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radiation pressure and the solar wind move the dust and ionized particles away from the nucleus, and a 
tail is formed. The composition of a comet is not well-known, but the nucleus is generally thought to be a 
“dirty snowball” made mostly of frozen water and organic and silicate compounds. When the comet is 
close enough to the Sun that the heat causes these volatile materials to boil off and carry solid particles 
with them, the comet is said to be “outgassing.” Consequently, force is exerted on the comet, and the 
resulting perturbations to the orbit should be taken into account when attempting to calculate a trajectory 
precisely enough to determine if a collision with Earth will occur. 

The most widely accepted method for modeling force due to outgassing was developed by Marsden, 
Sekanina, and Yeomans, as described in reference 7, and employed in reference 8, based on the assump- 
tion that the comet is an icy conglomerate — an object consisting mostly of frozen water that holds 
together bits and pieces of rock-like material. The contribution of outgassing to force per unit mass is 
expressed as 


f g - g(f)[-Aiai + A 2 a 2 + A 3 a 3 ] 
where the function g(r) of heliocentric distance r is given by 


g(r) = a 
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The parameters other than r are constants obtained from studies of vaporization rates of comet nucleus 
material. The normalizing constant a is defined such that g(l) = 1, and vq is the scale heliocentric 
distance of high outgassing activity with a value equal to 2.808 au for frozen water. The remaining 
constants are reported to be m = 2.15, n = 5.093, and l = 4.6142; substituting these values into equa- 
tion (65) and solving for the normalizing constant yields a = 0. 1 1 13. 

A right-handed, mutually perpendicular set of unit vectors a 1? a 2 , and a 3 is defined such that a 1 is in 
the direction of the position vector from S to C, a 2 lies in the orbital plane, and a 3 is normal to the orbital 
plane in the direction of the specific inertial angular momentum of C relative to S. Numerical values of 
the constant coefficients A ( and A 2 are calculated by studying changes in comet orbital periods, and are 
listed for more than 20 comets in table I of reference 7. The coefficient A 3 is usually neglected because 
the force normal to the orbit plane has no detectable effect on the period, therefore values for A 3 are not 
given. However, a normal force will influence a comet’s final Earth-encounter distance, so this 
component must be included in simulations involving prediction of a collision. 

Before undertaking such simulations, it is important to understand the magnitude of the perturbation to 
the comet’s trajectory caused by outgassing. This is accomplished with the aid of dynamical equations 
formed by adding the perturbing force per unit mass to the right hand member of the equations governing 
two-body motion, 


d~ n 

— 3 -r=-^-r + g(r)[A 1 a 1 + A 2 a 2 + A 3 a 3 ] (66) 

dt r 

where r is the position vector from S to C, and N d 2 r/dt 2 denotes the second derivative in frame N of r 
with respect to time. By comparing numerical solutions obtained through integration of equations (66) 
with [g(r) ^0] and without [g(r) = 0] the contribution of the force of expelled gas, one can quantify its 
effect on the time history of r. Integration is performed in connection with the hypothetical LPC 
whose orbital elements are given at the conclusion of the material dealing with range and range- 
rate measurements. Equations (66) are integrated from t Q to t,, at which times r = 5 au and 1 au, 

o 

respectively. The maximum values of the coefficients reported in reference 7 are used: A 1 = 3.61 au/10 
day 2 = 5.4 km/day 2 , and A 2 = 0.3269 au/10 8 day 2 = 0.49 km/day 2 . Because no values for A 3 are given, 
A 3 is set equal to A 2 . 

Figure 22 includes time histories of the differences in r-s ( - (z = 1, 2, 3) caused by outgassing, as well 
as the difference in the magnitude of r, which is approximately 4730 km by the time r = 1 au. However, 
the more important result is that differences in position are negligible until about 250 days, when r is 
approximately 2.2 au; therefore, a two-body trajectory furnishes a good reference orbit until this point, 
warranting the linearization about such a reference orbit performed in constructing the state transition 
matrix of equation (34). The state transition matrix should account for outgassing if the filter is to be used 
with observations of a comet made with r less than, say, 2.5 au, but need not do so for r greater than this 
number. 

The uncertainty that results from outgassing is accounted for with a white state noise of zero mean and 
time-varying diagonal covariance [(2(t/)] = [Qj], the elements of which are determined by assuming the 
perturbing force per unit mass is constant during the interval between observations. The uncertainty in 
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Figure 22. Change in position due to outgassing for a hypothetical comet. 


the force per unit mass is reflected in the standard deviations of the coefficients listed in reference 7: 
Oj = 2.096 km/day 2 , and a A = o A = 0.085 km/day 2 . The elements of [ Qj ] are given by 


QrS t j) = 


t \2 


n2 


(* = 1, 2, 3) 


(67) 


{[fg(^)-S(_3](ty tj- 1)} 


O' = 4, 5, 6) 


( 68 ) 


where f g (tj) is evaluated at tj with o A . in place of ( i = 1, 2, 3) in equation (64). As heliocentric 
distance r decreases, the elements of [Qj] grow along with g(r), as do the elements of the Kalman gain 
matrix [Kj], and the filter places increasing emphasis on the measurements in determining the best 
estimate of the state deviation. 

Once a best estimate { X'j } of the orbital elements is obtained from the sequential filter, an erroneous 
predicted miss distance is calculated by using the elements of (A' } as initial values for the differential 

equations (66) and by performing numerical integration with the MATLAB® function ODE 4 5 from the 
time tj of the final observation to the designed time of collision t k . 
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As it is not possible to know ahead of time the values of the outgassing coefficients for a newly 
discovered comet, a Monte Carlo analysis is performed with a set of random values uniformly distributed 
between the minimum and maximum values of A l and A 2 reported in reference 7. Random values for A 2 
are obtained from the same distribution used for A 2 . The results of this analysis are presented next. 

The Body-Targeting Plane 

Body-plane (B-plane) targeting is a method commonly used in the preliminary design of interplanetary 
missions to determine a distance of closest approach of a spacecraft to a celestial body of interest. A miss 
distance B calculated with this method can furnish a metric for quantifying orbit determination accuracy 
in connection with hypothetical comets that collide with Earth. The uncertainty in B can be expressed in 
terms of error ellipses constructed on the B-plane and is actually the more important result when 
analyzing the accuracy with which an orbit has been determined. 

A vector B, and a plane containing B, are constructed as follows. First, one neglects the gravitational 
force exerted by the Earth E on the hypothetical object C; therefore, E does not perturb the heliocentric 
orbit of C, and the velocity v in N of C relative to the mass center E* of E is regarded as constant in 

F IF* 

N during the time C is within the Earth’s sphere of influence. The magnitude of v is denoted by Voo. 

A unit vector S having the same direction as v is shown in figure 23, and is often considered to be 
parallel to the incoming asymptote of the geocentric hyperbolic trajectory (assuming C is not captured) 
resulting from gravitational attraction of the Earth that must in fact exist. In the absence of gravitational 
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force exerted by E, C continues along the incoming asymptote. The unit vector S can be obtained with 
the expression 


„C7E* 

S = ^ 


y c/s \t k )-y E * /s *(t k ) 


(69) 


r' /o* 

where v (f fc ) and v (fy.) are, respectively, the velocities in N of C and E*, relative to the Sun’s 
mass center S*, at the designed time of collision t k . 


The B-plane is defined such that it contains the point E* and is perpendicular to S; the plane contains 
two unit vectors T and R that are mutually perpendicular to each other and to S. The unit vector T is 
chosen to be in the ecliptic plane; therefore, it can be obtained from the relationship 


T = 


Sxs^ 


(Sxs 3 )-(sxs 3 ) 


n 1/2 


(70) 


where, as one will recall, s 3 is perpendicular to the ecliptic plane. Unit vector R is then simply 


R = SxT (71) 

The vector B denotes the position from E* to the point at which the incoming asymptote intersects the 
B-plane, which is the point of closest approach when one neglects the gravitational attraction of E (or, 
equivalently, neglects the mass of E). The magnitude B ofB is the corresponding distance of closest 
approach. 

As previously mentioned, the geocentric velocity v = vjS is assumed to remain constant in N for 
the purpose of constructing B; therefore, C must necessarily travel on a straight line between the point of 
closest approach and the position at the designed time of collision, as shown in figure 24. Hence, the 
specific angular momentum h in N of C relative to E* is given by 


h = B x v JS = 8 x v JS 


(72) 


where 8 is the position vector from E* to C at the designed time of collision, obtained from the process of 
orbit determination. From this relationship, an expression for B can be derived through premultiplication 
with S in a cross product, 


Sx h = Sx (e x vjS) = Sx (B x vjS) 

= v„[(S-S)B-(SB)S] 

= v B 


because S - S= 1, and S - B = 0 by definition. Hence, 


B= S x h = sx(exS) 


(73) 


(74) 
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Figure 24. Erroneous predicted miss position, and closest approach. 


and, because B lies in the plane formed by R and T, it can always be expressed as 

B = B R R + B r T (75) 

As described previously, the orbit of C is designed to pass through the ecliptic at a heliocentric 
distance of 1 au, a point that is assumed to be coincident with E*. A perfect determination of the orbit 
must necessarily yield e = B = 0; however, inevitable random noise present in observations will lead to a 
prediction that e and B each differ from 0. If the magnitude of either of these position vectors is 
determined to be, say, less than or equal to 1 Earth radius (R£), one would be justified in predicting a 
collision unless there is a large uncertainty in the orbit solution as quantified by the covariance of e or B, 
which bears the following relationship to the estimate error covariance matrix [Pj\ in equation (31). 

A time update according to equation (26) is performed in order to obtain the covariance matrix at the 
designed time of collision t k , from the covariance matrix \Pj\ after a measurement update at tj, 

[p k ]=mt k ,tj)[Pj]mt k ,tj)] T +[Qj ] ( 76 ) 


The effects of outgassing can be reflected in [QA as set forth in equations (67) and (68) for short 
intervals of time. Alternatively, one may divide a large time interval into shorter ones over which f„ 
is regarded as constant in N. When 1 year of warning time separates t k and tj , dividing the interval into 
20 equal parts ensures that [Qy] is not dependent on the propagation time step and that the constant 
outgassing force assumption is still valid. The covariance time update for the shorter intervals (each 
approximately 2 weeks in length) is modified, 

[p j+1 ]=mt j+l ,tj)][p j ]mtj +1 ,t j )f +[Qj i (77) 


to employ the average value of the state noise covariance 

IQj] i 1Gjl + ‘ Qj+l1 


(78) 
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Equation (77) is applied successively until fy +1 becomes the time of collision t k , yielding a 6 x 6 
covariance matrix [P k ]. 

The upper left 3x3 partition of f 7). ] deals with position and is therefore of the most interest in what 
follows; it is denoted by the symbol [F r ]. As indicated in connection with the modeling of comet 
outgassing, the state transition matrix should be constructed according to reference 8 when outgassing 
becomes a significant perturbation; in lieu of this, [F r ] can be supplemented by adding a covariance 
matrix assembled from a difference in position associated with uncertainty resulting from outgassing. 
Equation (66) is numerically integrated from the time of the first observation t\ to the time of collision t/ ( , 
first with fg = 0 , and then three times, with one of a A used in place of A f (i = 1, 2, 3) in equation (64). 
The difference in r (t k ) obtained with c A ., and with f g = 0 , is expressed in terms of unit vectors S] , §2, 
and S3 and presented as a column matrix {8} ( - , allowing one to form the 3x3 covariance matrix 

3 

[^] = [^]+I(5} ; {6}f (79) 

i= 1 

This matrix (like the rest of [/),]) is associated with unit vectors S3, $2, and S3, but it is more convenient 
to work with a covariance that describes the uncertainties in the directions marked by S, T, and R. The 
transformed covariance, indicated with [Fg], is formed as 

[F g ] = [ 5 C B ] T [F 5 ][ 5 C fl ] (80) 

where the elements of the direction cosine matrix [ C ] are defined as C,q = s -S, C i2 = s, T, and 
S cf 3 =Si- R (i = l,2, 3). 

As discussed in section 4.16 of reference 5, the function 

{{x k }-{x k }f[P B r\{x k }-{x k })=l 2 (81) 

describes a three-dimensional ellipsoid, where the portion of the state deviation associated with position 
at t k is {x k }, having covariance [Fg] and mean {T fc }, and where / is a constant. The ellipse correspond- 
ing to a given value of /, formed by the intersection of the ellipsoid with the B-plane (shown in fig. 23), is 
obtained from the 2 x 2 partition of [Fg] associated with T and R, represented by [/),], which can be 
diagonalized 


[P p ] = [ B C F ] T [ P b ][ B C P ] 


(82) 


BP — 

where [ C ] is a matrix containing the normalized eigenvectors of [F^], with elements defined as 

B PA" ~ B P A ^ ^ — 

Cu = T p, and °C r 2i = R p, (i = l, 2, 3). The diagonal covariance matrix [ P p ] contains the eigen- 
values of [P/.,] associated with the eigenvectors of unit length pi and p 2 , which are parallel to the 
principal axes of the error ellipse in the B-plane. The orientation of the principal axes of the ellipse with 
respect to T and R can be obtained as 0 = arctan C21/ Cn). The square roots of the two nonzero 
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elements of [ P p ] are the semimajor and semiminor axes of the ellipse, denoted as Oi and g 2 , respectively. 
The difference in time that it takes for the comet to travel between the positions marked by B and e is 
assumed to be negligible, and the error ellipse therefore describes the uncertainty in the estimate that the 
comet will pass through the position indicated by B. The bivariate Gaussian probability density function 
can be expressed in invariant form according to equation (2.2-39) of reference 9, or equation (A.19.1) of 
reference 5, 


fix i,x 2 ) 


1 -qOM-WlVr 1 ({*}-{*}) 


2n\P\ 


1/2 


(83) 


where [i 5 ] can be either of the covariance matrices [TJJ or [ P p ], \P\ denotes the determinant of [P], and 

r /v /v 

{x} = [xj x 2 \ contains the random variables of the state deviation at t k associated with T and R in 
the case of [TJ,], or with pj and p 2 i n the case °f [-Pp]- The random variables have a mean of {x}, a 

column matrix containing the two dot products B- T = Bj and B R = B R , or B p | and B p 2 . Contours 
of constant / are ellipses. 

Inspection of results from a Monte Carlo analysis shows the validity of the assumption of linearity 
inherent in the development of the sequential fdter. The object of interest is the hypothetical comet 
described previously at the conclusion of the section dealing with range and range-rate measurements, 
without the perturbation induced by outgassing. The extended Kalman fdter is employed to produce 1 00 
orbit solutions from 91 optical observations spaced 1 day apart, made from a single observatory whose 
resolution p is 0.01 arcsec. Each solution begins with a preliminary orbit determined by the method of 

Gauss, and processes a different set of measurements of longitude and latitude, whose random errors are 
formed according to step 1 of the procedure set forth for obtaining a priori information. The procedure 
yields 100 values of B and associated covariance matrices | P/,} (or [ P p ]), each of which corresponds to 
an ellipse whose center is the position given by B and indicated with a point in figure 25. The average 
values of 01 , o 2 , 0, Bp, and Br are given in the second column of table 2. Counterparts to the first three 

of these parameters may be obtained by forming the covariance of the 100 state deviations {x} produced 
by the filter after the final measurement update, and mapping this single covariance matrix to t k . Error 

ellipse parameters obtained in this way are given in the third column of table 2, and the la ellipse 
centered on the average values of B p and Br is shown in figure 25 to be a reasonable representation of the 
distribution of the individual intersection points. 


The time of closest approach t h can be determined from the designed time of collision q. by referring 

to figure 24 and recalling the assumption that the geocentric velocity v = v m S is assumed to remain 
constant in N for the purpose of constructing B. The time tp — tp required for C to travel from the 
position marked by B to the position marked by e is given by 

t k -t b = —(e- B)-S=— (84) 

V V 

' 00 T 00 

The standard deviation c t of this time interval obtained from the 100 orbit solutions is given in the fourth 
column of table 2. 
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Figure 25. la error ellipse for long-period comet. 
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Figure 26. Variation in 1 a error ellipse for long-period comet. 
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Alternatively, the accuracy of the estimate of the time at which the comet will arrive at the position 
marked by B is indicated by the extent of the three-dimensional error ellipsoid normal to the B-plane. 
Letting denote the square root of the element of [P B ] associated with S, the uncertainty in time of 
arrival at the B-plane is given simply by 


0 / = c 5 /v oo (85) 

Table 2 contains the average values of 05 and time of arrival standard deviation in the second column 
and the values obtained from the single covariance matrix in the third column. Values of 0 | , 0 2 , 05 , or 
a T corresponding to 20 and 3o ellipses are calculated easily by multiplying the lo value in table 2 by 
2 or 3. 

The similarity of 0 f in the fourth column to the values recorded in the second and third columns 
shows the errors in time of arrival have the same standard deviation as the travel times between the 
positions marked by B and 8. The close agreement in the values of all parameters in the second and third 
columns of table 2 is one indication the assumption of linearity made in deriving the sequential filter is 
justified. Additional support is furnished in figure 26: the individual error ellipse parameters for each of 
the 100 Monte Carlo simulations are virtually the same, with only slight variations distributed on either 
side of the mean (shown with a solid blue line). The similarity in these error ellipse parameters is 
evidence of linearity; in other words, the covariance of the state deviation at the designed time of collision 
is virtually independent of the point at which the orbit intersects the B-plane. Consequently, there is no 
need to perform additional Monte Carlo analysis. 

The foregoing conclusions remain unchanged after including the effect of outgassing, as indicated in 
table 3. A comparison of tables 2 and 3 reveals that comet outgassing increases the size of the lo error 
ellipse, signaling a less accurate determination of the orbit. 


Table 2. lo Error Ellipse Parameters for Sample Comet 


Parameter 

Average 

Single 

Eq. (84) 

0b Re 

0.072535 

0.071114 


02> Re 

0.038867 

0.045221 


9, deg 

-16.579420 

-8.787936 


B t , Re 

0.008113 



Br , Re 

0.002959 



0 * r e 

0.078602 

0.077102 


a t , min 

0.231822 

0.227399 

0.227400 

Table 3. 

la Error Ellipse Parameters for Sample Outgassing Comet 

Parameter 

Average 

Single 

Eq. (84) 

0b Re 

0.287882 

0.287184 


02> Re 

0.064515 

0.069989 


9, deg 

-31.912291 

-31.814645 


Bt, Re 

-0.210096 



Br, Re 

-0.108374 



0 * r e 

0.334448 

0.334101 


a t , min 

0.986446 

0.985425 

0.774256 
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Probability Associated With B 

The probability that an object passes through the B-plane within a certain distance of the position 
indicated by B is found by integrating the probability density function over the area of a circle. The result 
can be used in two ways, provided the assumption of linearity is warranted. First, in the unique 
circumstance that occurs because the orbit of our hypothetical object is known precisely, the probability 
quantifies the accuracy of the orbit determined by the sequential filter. Second, for the general case in 
which the true orbit is unknown, a probability of collision between the object and Earth can be obtained. 

In what follows, numerical examples are created with the hypothetical comet described in the 
discussion regarding range and range-rate measurements; again, the filter processes 91 optical observa- 
tions from a single observatory, taken 1 day apart, and measurement errors are normally distributed, based 
upon a resolution of p = 0.01 arcsec. In view of the demonstration in the preceding section, the 
assumption of linearity is justified, and we dispense with Monte Carlo analysis; furthermore, the 
eigenvalues and eigenvectors of the covariance matrix [P b ] are considered independent of B (the mean of 
the state deviation) obtained from the sequential filter in the orbit determination process. The actual 
orbital parameters of the hypothetical object are adopted to produce the reference trajectory for the 
sequential algorithm, and a conventional filter is put into service so that the reference orbit remains 
unaltered. 


The bivariate probability density function given by equation (83) is integrated numerically over the 
area of a circle whose center is given by the coordinate pair (T] ,x 2 ), the elements of which are used to 
form the matrix { x}. The center of the circle varies over several points in a coarsely spaced grid in the 
neighborhood of (0, 0) in the B-plane. For each center point, the function /(jq,jt 2 ) is evaluated at every 
point of a fine, evenly spaced grid that approximates a circle, and then multiplied by the area enclosed by 
four neighboring grid points (arranged in a square); the products thus obtained are summed to yield the 
probability that the object passes through the B-plane, somewhere within the circle. 


The radius of the circle is not critical for the purpose of gauging orbit determination accuracy. 
However, a particular radius is required for calculating probability of collision; the area of the circle must 
be equal to the effective collision cross section of the planet, which is larger than the physical cross 
section in order to account for Earth’s gravitational attraction that has been ignored up to now. The 
effective collision radius r c is given by equation (8.3-30) in reference 1, 




2 I X E 

Re 


( 86 ) 


where R E is the physical radius of Earth (regarded as a sphere), and \i E is the Earth’s gravitational 
parameter. 

The probability that the object passes through a circle of radius r c centered at the Earth’s mass center 
E* , (0, 0), is found to be 99.9 percent. The probability of passing through circles with other centers is 
likewise computed, and each concentric ellipse in figure 27 is a locus of points for which the probability 
is equal to 1 in 10 z (z = 1, 3, 6). Because the actual orbit is designed to result in B = 0, it can be said that 
the higher the probability, the more accurate the orbit determination. With the luxury of knowing the 
actual orbit corresponds to e = B = 0, the information presented in figure 27 gives reason to be 100 times 
more confident in a filter solution when B marks a point on the contour labeled -1 than when it marks a 
point on the -3 contour. 
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Figure 27. Probability of passing through a circle of radius r c . 


In addition to gauging orbit determination accuracy, figure 27 can be used to indicate probability 
of collision. In the event the filter solution yields B = 0 (rather unlikely), the probability of collision is 
99.9 percent, and the probability is one in a million that the object crosses the B-plane through a circle of 
radius r c centered at, say, (1, -1). On the other hand, if the filter solution gives B = 1T-1R (a more 
likely result than B = 0), the probability that the comet passes within a distance r c of this point is 
99.9 percent; furthermore, there is one chance in a million that the comet will instead cross the B-plane at 
(0, 0) and collide with Earth. This interpretation of figure 27 rests on the assumption of linearity and the 
concomitant concept that the uncertainty in B is described by the same error ellipse, no matter what the 
direction and magnitude of B. In other words, the same error ellipse can be translated or shifted to any 
point on the B-plane. 

The presence of outgassing can be expected to reduce the probability of passing through B = 0; 
measurements of range and range-rate can be expected to increase it. The orbit solution based on angular 
measurements with p = 0.01 arcsec from 91 observations at a single observatory is so accurate that neither 
of these effects are observed; however, varying the number of observations brings them to light in fig- 
ure 28. The plot at the top shows the increase in central probability with number of observations. The 
second and third plots show the respective collision probabilities associated with solutions of B = IT and 
B = 1R. There are six curves in each plot corresponding to three combinations of observation types. The 
first consists of two angular measurements at each observation, the second includes range measurements 

'J 

at a resolution of 10 km, and the third adds measurements of range-rate at a resolution of 1 m/s. For 
each combination, results with outgassing are compared with those in which outgassing is absent. 
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Center of Earth 




Figure 28. Probability, outgassing, and measurements of range and range-rate. 


Batch Filter for Near-Earth Asteroids 

Having examined orbit determination of LPCs in some detail, we turn our attention now to potentially 
dangerous objects whose orbital parameters are rather different, NEAs that are so named when they have 
perihelial distances less than 1.3 au. These asteroids are further classified by dividing them into groups 
based upon values of semimajor axis a, aphelial distance r a , and perihelial distance r p . Earth-crossing 
asteroids with a < 1 au and r a > 0.983 au are members of the Aten class. Earth-crossing asteroids with 
a > 1 au and r p < 1.017 au belong to the Apollo class. Earth-approaching asteroids with a > 1 au and 
1.017 < r p < 1.3 au are classified as Amors (these asteroid orbits do not actually cross Earth’s orbit but 
come close enough that a perturbation could cause a collision). It is possible that there exists another 
class of Earth-approaching asteroids for which r a < 1 au, although no asteroids of this type have been 
detected to date. 

The following discussion of analysis is facilitated by the definition of two types of Earth-crossing 
asteroids that are more general than the foregoing classes. “Interior” asteroids are said to be those with 
r p < 1 au and r a = 1 au, while “exterior” asteroids have r p <1 au and r a > 1 au. Figures 29 and 30 
contain examples of these two types of orbits, where the Earth E is shown in a circular orbit of radius 1 au 


71 



Figure 29. Example of interior asteroid orbit. 




Figure 30. Examples of exterior asteroid orbits. 


72 


about the Sun S, and an asteroid A is shown in a coplanar elliptical orbit. The negative signs of eccentric 
anomalies E and E k required for preperihelial collision with LPCs (see the discussion of warning time) 
apply also to exterior asteroids; however, collisions of interior asteroids must always take place after 
perihelion, therefore the sign changes are unnecessary. 

The focus of the present analysis is upon asteroids discovered less than 1 orbital period prior to 
collision. Preliminary orbit determination for interior and exterior asteroids is performed in much the 
same way as it is in connection with LPCs, with the following details regarding the seven-step procedure 
given in connection with erroneous predicted miss distance. 

In the first step, the earliest opportunity for observations t 0 is specified as 60 days prior to the time of 
collision t k , and the initial value of true anomaly v(?q) is calculated accordingly. The reason for basing 
v(fo) on time to collision rather than a specified value of r has to do with the observation times chosen in 
the third step. Four observations permit the determination of orbits inclined 0° or 180° to the ecliptic; we 
wish to make the fourth and final observation at least 2 weeks prior to collision in order to provide some 
amount of warning. Hence, the times q, t 2 , t 2 , and f 4 of the observations are taken to be 11, 22, 33, 
and 44 days after t 0 , leaving 1 6 days between the final observation and collision. A single observatory in 
a circular heliocentric orbit of radius 1 au in the ecliptic plane is used to make the observations, and the 
initial true longitude L 0 takes on values of 0°, 90°, 180°, and 270°. 

The uniformly distributed measurement errors associated with the fourth step are bounded in 
magnitude by a telescope resolution of p = 0. 1 arcsec, and once again are constructed with a pseudo- 
random number algorithm. Average values e of erroneous predicted miss distances for 100 sets of 
measurements are shown. 

In connection with the fifth step, the position and velocity at q are to be determined by the method of 
least squares. It so happens that the algorithm is quite sensitive to the initial estimate {T} of position and 
velocity; therefore, values corresponding to day 11 are used for exterior asteroids, and to day 12 for 
interior asteroids. The MATLAB® algorithm failed to converge for a number of measurement sets in a 
handful of cases; nevertheless, a value of £ is computed from the unconverged estimates. The sensitivity 
to the initial guess, and the lack of convergence, may be due to changes in an asteroid position that are 
significantly larger than changes in a comet orbit over the same interval of time. 

Orbit determination for interior asteroids is discussed first, and exterior asteroids follow. 

Interior Asteroid Orbits 

Classical orbital elements of interior asteroids are constructed with r a = 1.0 au, r p = 0.2, 0.4, 0.6, and 
0.8 au, i = -40°, -30°, -20°, ..., 40°, and Q = 0°, 45°, 90°, 135°, and 180°. For each of these 180 sets of 
orbital elements, the argument of perihelion co is determined to be 180°, according to equation (5). The 
initial true longitude L 0 of the observatory takes on four values, each differing by 90°, resulting in a total 
of 4 x 180 = 720 cases. 

The results from this analysis are presented in figure 3 1 as a function of asteroid orbital elements i and 
r p \ each data point represents £ for one case, and a family of 20 cases (in which Q and L 0 vary) form 
one stack of points. At the lower values of r p one can see an inverse relationship of £ to i, which is not 
unexpected because orbits that are coplanar with the observatory require a minimum of four observations 
to be determined. Although not as evident at this scale, the relationship exists at the higher values of r p 
as well. At i = 0° and r p = 0.6 au one notices two outlying values of £, approximately 0.45 lunar 
distance. These two cases are distinguished from each other only by the value of L 0 (and of Q, which is 
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Figure 3 1 . Interior near-Earth asteroids. 


undefined for i = 0). For a third case in this family the least squares algorithm was unable to converge on 
a solution for 38 of the 100 measurement sets, and the associated data point is not shown. Of the 719 
results displayed, 697 values of e are less than 0.05 lunar distance, or about 3 Earth radii. When one 
holds i constant and plots the results in two dimensions, e is seen to vary proportionally to r p when 
i = ±40° and inversely to r p when i = ±10°. When displayed as functions of orbital elements other than 
i and r p , no remarkable trends are visible in the results. 

In eight other cases a certain number of measurement sets did not produce convergence. These cases 
have common values of r a = 1 au, r p = 0.2 au, and Q - L 0 - 45°. The nine cases that suffered from 
convergence problems are listed in table 4, together with their orbital parameters, initial true longitude of 
the observatory, and the number of trials out of 1 00 that did not converge. 


Table 4. Interior Asteroid Orbits With Unconverged Trials 


Case 

r a , au 

r p , au 

i, deg 

Q, deg 

to, deg 

Vo, deg 

To, deg 

Trials 

115 

1.0 

0.6 

0 

180 

180 

123.72 

180 

38 

362 

1.0 

0.2 

-40 

45 

180 

132.52 

0 

85 

367 

1.0 

0.2 

-30 

45 

180 

132.52 

0 

3 

397 

1.0 

0.2 

30 

45 

180 

132.52 

0 

4 

402 

1.0 

0.2 

40 

45 

180 

132.52 

0 

86 

544 

1.0 

0.2 

-40 

135 

180 

132.52 

90 

1 

549 

1.0 

0.2 

-30 

135 

180 

132.52 

90 

10 

579 

1.0 

0.2 

30 

135 

180 

132.52 

90 

9 

584 

1.0 

0.2 

40 

135 

180 

132.52 

90 

7 
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Average e, lunar distance Average £, lunar distance 




Figure 33. Accuracy of exterior NEA orbits, as a function of i and r p . 
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Exterior Asteroid Orbits 


For exterior asteroid orbits, 1350 sets of classical orbital elements are constructed with r a = 1.5, 3.0, 
4.5, ..., 9.0 au, r p = 0.2, 0.4, 0.6, .... 1.0 au, i = -40°, -30°, -20°, ..., 40°, and Q = 0°, 45°, 90°, 135°, 
and 180°. For each set, the argument of perihelion co is computed according to equation (5). L 0 takes on 
the same four values as in the interior cases; thus, a total of 4 x 1350 = 5400 cases for exterior asteroid 
orbits are examined. 

The average values £ of erroneous predicted miss distance are shown in figures 32 and 33 as functions 
of i, and r a and r p , respectively. Relatively poor orbit determination is exhibited once again for asteroids 
traveling in orbits that have low inclination with respect to the observatory, when only four observations 
are obtained. Accuracy generally decreases for smaller values of r p and for larger values of r a . 

Of the 5400 cases examined, 5389 yielded £ less than 1 lunar distance. In two of the remaining cases 
£ was reported to be approximately 4.8 lunar distances; however, the least squares algorithm was unable 
to converge on a solution for 15 and 38 of the 100 measurement sets, and the associated data points are 
not shown. The largest value of £ shown in figures 32 and 33 is 2.14 lunar distances, and the smallest is 
0.00125 lunar distance, or 0.08 of an Earth radius. In seven other cases involving an object whose orbit is 
coplanar with that of the observatory (/ = 0), a certain number of measurement sets did not produce 
convergence, as reported in table 5. 


Table 5. Exterior Asteroid Orbits With Unconverged Trials 


Case 

r a , au 

r p , au 

i, deg 

Q, deg 

co, deg 

Vo, deg 

To, deg 

Trials 

564 

4.5 

0.6 

0 

135 

85.59 

-120.36 

180 

15 

2991 

3.0 

0.4 

0 

0 

112.62 

-140.51 

0 

49 

3082 

3.0 

0.8 

0 

45 

62.96 

-106.90 

0 

5 

3307 

4.5 

0.8 

0 

45 

59.10 

-102.53 

0 

32 

3532 

6.0 

0.8 

0 

45 

57.42 

-100.64 

0 

30 

3757 

7.5 

0.8 

0 

45 

56.48 

-99.58 

0 

24 

3982 

9.0 

0.8 

0 

45 

55.88 

-98.91 

0 

19 

4479 

3.0 

1.0 

0 

135 

0 

-62.90 

0 

1 

4612 

4.5 

0.6 

0 

45 

85.59 

-120.37 

0 

38 


Although the analysis for NEAs has not been extended to include multiple observatories or improved 
telescope resolution, it is likely that configurations of two or more observatories, together with multiple 
observations, will improve preliminary orbit determination for asteroids, as has been demonstrated in 
connection with LPCs. 

Concluding Remarks 

An examination has been made of the effects of several factors upon determination of orbits of comets 
and asteroids on a collision course with Earth, including the time interval between successive observa- 
tions, the distance at which the object is detected, the resolution of the telescope, the orbital parameters of 
the object, the number of observations, and, in a rather limited manner, placement in the solar system of 
one or more observatories. The primary factor is seen to be the length of the data arc. As mentioned at 
the outset, the analysis is based on a crucial assumption that optical measurements can in fact be obtained 
from a telescope; the validity of this assumption depends upon several factors that have not yet been 
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addressed, among them the quality of the instrument, the extent of solar illumination of the object (solar 
phase angle), the optical reflectivity (albedo) of the object, occultation by the Sun or other bodies, and 
interference of light from sources such as the stars in the galactic plane, zodiacal dust, and the Sun. 

Provided the measurements are available, the following conclusions are made based upon a study of 
many hypothetical objects, each of which has been given an orbit resulting in a collision with Earth at a 
specified time and a specified position. As used here, accuracy is expressed in terms of erroneous 
predicted miss distance, the difference between an object’s determined position and the specified position, 
evaluated at the specified time of collision. A single observatory in a circular heliocentric orbit of radius 
1 au in the ecliptic plane, possessing an optical resolution of 0.1 arcsec, provides good preliminary orbit 
determination accurate to less than 8 lunar distances for long-period comets (LPCs) observed over a 
66-day arc at a heliocentric distance of 6.5 au. Multiple observations (approximately 100, spaced 1 day 
apart) improve the accuracy to less than 14 Earth radii, and a second observatory phased 180° from the 
first improves the accuracy further still to less than 1 Earth radius. Of the multiobservatory config- 
urations examined here, two observatories of equal resolution in circular heliocentric orbits of radius 1 au 
in the ecliptic plane, phased 180° apart, provide the best balance of cost and accuracy in orbit 
determination for LPCs. Thirty to forty observations taken over a 98-day period appear to give results 
nearly as good as 99 observations taken 1 day apart. For near-Earth asteroids observed over an interval 
between 60 and 16 days before collision, the accuracy is better than 2 lunar distances for the exterior 
variety and 0.5 lunar distance for interior asteroids. 

The results reported here give a preliminary indication of what is required to predict impending 
collisions accurately. We must keep in mind that such predictions must be reliable enough to support a 
decision to expend considerable resources in attempting to alter the orbit of a dangerous body; very 
accurate knowledge of the orbit is probably essential if orbital modification is to steer the object away 
from a collision, rather than toward one. Moreover, predictions must have a degree of integrity that 
virtually eliminates false alarms and unfounded panic. The results indicate that it may be possible to 
make reliable forecasts with two observatories whose angular resolution is on the order of 0.1 arcsec, or 
with a single observatory whose resolution is better by 1 to 3 orders of magnitude. Additional study of 
these two alternatives will have to weigh the advantages in geometry and redundancy of multiple 
observatories against the expense of putting them in place and maintaining them; a single observatory 
near Earth could be easier to maintain, but less likely to be in the best position for obtaining 
measurements. 

Outgassing does not perturb a comet’s orbit appreciably until it reaches a heliocentric distance of 
approximately 2.5 au; the perturbation in position is less than the Earth’s radius by the time the comet 
reaches 1 au. Analysis of dangerous LPCs performed with Kalman filters reveals several interesting 
results. With a mix of weekly and daily optical observations from one or two observatories, a credible 
warning of collision, in the absence of outgassing, can be expected at least a year in advance if the object 
is spotted when it is about 6.5 au from the Sun. If a comet is discovered at a heliocentric distance of 5 au, 
then supplementing daily angular measurements whose resolution is 0.01 arcsec, with range measure- 
ments accurate to 1 000 km and range-rate measurements good to 1 m/s, is beneficial when dealing with a 
short data arc of 15 to 30 days, but doing so offers little advantage when the arc is 80 or 90 days long. 
The assumption of linearity upon which the sequential filters depend is warranted, even in the presence of 
outgassing, when daily angular measurements with a resolution of 0.01 arcsec are available for 3 months. 
These sorts of measurements yield extremely accurate orbit determination, even with only 2 months of 
observations; the addition of range data can reduce the necessary data arc to 1.5 months, and only 
1 month is needed when range-rate is included. 
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Future work should investigate the effects of perturbations that have been ignored here, including the 
gravitational attraction of Jupiter, Saturn, and other planets, solar radiation pressure, and relativistic 
effects. The advantages of placing observatories in heliocentric orbits other than those studied here 
deserves careful consideration. In addition, the consequences of missed observations (due to the reasons 
stated earlier) should be examined. 
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